FE_Publication_Eukaryota_Biomass
library(ROBITaxonomy)##
## Attachement du package : 'ROBITaxonomy'## L'objet suivant est masqué depuis 'package:stats':
##
## familylibrary(ROBITools)## Warning: remplacement de l'importation précédente 'dbplyr::ident' par
## 'dplyr::ident' lors du chargement de 'ROBITools'## Warning: remplacement de l'importation précédente 'dbplyr::sql' par 'dplyr::sql'
## lors du chargement de 'ROBITools'## Warning: remplacement de l'importation précédente 'dplyr::union' par
## 'igraph::union' lors du chargement de 'ROBITools'## Warning: remplacement de l'importation précédente 'dplyr::as_data_frame' par
## 'igraph::as_data_frame' lors du chargement de 'ROBITools'## Warning: remplacement de l'importation précédente 'dplyr::groups' par
## 'igraph::groups' lors du chargement de 'ROBITools'## Warning: remplacement de l'importation précédente 'ROBITaxonomy::path' par
## 'igraph::path' lors du chargement de 'ROBITools'## Warning: remplacement de l'importation précédente 'ROBITaxonomy::parent' par
## 'igraph::parent' lors du chargement de 'ROBITools'## ROBITools packagelibrary(vegan)## Le chargement a nécessité le package : permute## Le chargement a nécessité le package : lattice## This is vegan 2.5-7##
## Attachement du package : 'vegan'## L'objet suivant est masqué depuis 'package:ROBITools':
##
## rarefylibrary(ade4)
library(ggplot2)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ tibble 3.1.2 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1
## ✓ purrr 0.3.4## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(zoo)##
## Attachement du package : 'zoo'## Les objets suivants sont masqués depuis 'package:base':
##
## as.Date, as.Date.numericlibrary(viridisLite)
library(viridis)
library(ggpubr)
library(grid)
library(gridExtra) ##
## Attachement du package : 'gridExtra'## L'objet suivant est masqué depuis 'package:dplyr':
##
## combine# The palette with grey:
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
# The palette with black:
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")taxo = read.taxonomy("Data/ncbi20210212")The Read contingency table
reads = read.csv("Data/Faeces/FE.Euka.samples.reads.csv",
header = TRUE,
row.names = 1)
reads = as.matrix(reads)
reads = decostand(reads,method = "total")The sample description table
samples = read.csv("Data/Faeces/FE.Euka.samples.samples.csv",
header = TRUE,
row.names = 1)The MOTU description table
motus = read.csv("Data/Faeces/FE.Euka.samples.motus.csv",
header = TRUE,
row.names = 1)Create an object, where you merge the three tables
Euka01 = metabarcoding.data(reads = reads,
samples = samples,
motus = motus)motus.hist = colMeans(reads(Euka01))
Euka01@motus$mean_ref_freq = motus.hist
Euka01 = Euka01[,order(motus.hist,decreasing = TRUE)]motus(Euka01)[1:20,c("scientific_name","mean_ref_freq", "rank")]viridiplantae_taxid <- ecofind(taxo,"^Viridiplantae$")
Euka01@motus$is_plant <- is.subcladeof(taxo,Euka01@motus$taxid,parent = viridiplantae_taxid)| # Loading of metadata |
r library(readr) samples.metadata_biomass <- read_delim("Data/metadata.csv", ";", escape_double = FALSE, trim_ws = TRUE) |
## ## ── Column specification ──────────────────────────────────────────────────────── ## cols( ## sample_id = col_character(), ## Animal_ID = col_character(), ## Sample_number = col_double(), ## Date = col_character(), ## Sample_time = col_time(format = ""), ## times_from_birch = col_double(), ## Fed_biomass = col_double() ## ) |
r samples.metadata_biomass = left_join(Euka01@samples,samples.metadata_biomass,by="sample_id") %>% select(-name, -replicate) |
r Euka01@samples = samples.metadata_biomass rownames(Euka01@samples) = as.character(Euka01@samples$sample_id) |
| Timeshift to synchronize time 0 to 21.00, 23.january 2018 |
| ```r Euka01@samples[Euka01@samples$Animal_ID==“X”,“times_from_birch”] = Euka01@samples[Euka01@samples$Animal_ID==“X”,“times_from_birch”] + 6 |
| Euka01@samples[Euka01@samples$Animal_ID==“Y”,“times_from_birch”] = Euka01@samples[Euka01@samples$Animal_ID==“Y”,“times_from_birch”] + 3 |
| Euka01@samples[Euka01@samples$Animal_ID==“Z”,“times_from_birch”] = Euka01@samples[Euka01@samples$Animal_ID==“Z”,“times_from_birch”] + 4 ``` |
r Animal_id <- rep("Unknown", nrow(Euka01)) Animal_id[which(samples(Euka01)$animal_id == "X")] <- "9/10" Animal_id[which(samples(Euka01)$animal_id == "Y")] <- "10/10" Animal_id[which(samples(Euka01)$animal_id == "Z")] <- "12/10" Euka01@samples$animal_id <- factor(Animal_id, levels = c("9/10", "10/10", "12/10") ) |
r Euka01@samples$Fed_biomass = factor(as.character(Euka01@samples$Fed_biomass), levels=c("20","500","2000")) |
Animal = Euka01@samples$animal_idLecanoromycetes.taxid = ecofind(taxo,"^Lecanoromycetes$")
is.a.Lecanoromycetes = is.subcladeof(taxonomy = taxo,Euka01@motus$taxid,Lecanoromycetes.taxid)
is.a.Lecanoromycetes[is.na(is.a.Lecanoromycetes)]=FALSE
sum(is.a.Lecanoromycetes,na.rm = TRUE)## [1] 1Euka01@motus[is.a.Lecanoromycetes,c("scientific_name","mean_ref_freq")]Lecanoromycetes = Euka01@reads[,is.a.Lecanoromycetes]data_euka01 = cbind(as.data.frame(Euka01@reads),
time = samples(Euka01)$times_from_birch/24,
animal_id = samples(Euka01)$animal_id)
start_time=11
end_time=23
plot1 = ggplot(data = data_euka01,
aes(x = time,
y = EUKAP2_00000018,
col=animal_id)) +
geom_point() +
geom_smooth(data = data_euka01 %>% filter(time >= start_time & time <= end_time),
method = lm,show.legend = FALSE) +
scale_color_manual(values=cbbPalette) +
stat_summary_bin(fun = mean, geom = "line") +
scale_y_log10() +
geom_vline (xintercept = c(0.54,1.54,2.54), colour = cbbPalette[6]) +
geom_vline (xintercept = c(5.54,6.54,7.54), colour = cbbPalette[7]) +
geom_vline (xintercept = c(10.54,11.54,12.54), colour = cbbPalette[8]) +
geom_vline (xintercept = 26, colour = "black") +
guides(color=guide_legend(title="Individual")) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 15),limits = c(-2,30)) +
theme_minimal() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black")) +
labs(y="Proportion of Lecanoromycetidae DNA",
x="Time (days after Betula pubescens was fed)") +
labs(x=expression('Time (days after'~italic(Betula)~ ~italic( pubescens)~ 'was fed)'),
fill="Var1")
plot1
ggsave("Euka.Biomass_for_publication_oneline.png", plot = plot1, width = 25, height = 16, units = c("cm"))## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 8 rows containing non-finite values (stat_smooth).## Warning: Removed 37 rows containing non-finite values (stat_summary_bin).## Warning: Removed 15 rows containing missing values (geom_point).Estimation of half-life time
\[ \begin{array} t1 &= a \times log(x) + b \\ t2 &= a \times log(x/2) + b \end{array} \] half-life time is \(t2-t1\)
\[ \begin{array} t1 &= a \times log(x) + b \\ t2 &= a \times (log(x) - log(2)) + b \\ t2 &= a \; log(x) - a \; log(2) + b \\ \end{array} \]
\[ t2 - t1 = a \; log(x) - a \; log(2) + b - (a \times log(x) + b) \] \[ t2 - t1 = - a \; log(2) \]
data_euka01 %>%
filter(time >= start_time &
time <= end_time &
EUKAP2_00000018 > 0 &
animal_id=="9/10") %>%
lm(time ~ log(EUKAP2_00000018),data=.) -> lm_9_10
data_euka01 %>%
filter(time >= start_time &
time <= end_time &
EUKAP2_00000018 > 0 &
animal_id=="10/10") %>%
lm(time ~ log(EUKAP2_00000018),data=.) -> lm_10_10
data_euka01 %>%
filter(time >= start_time &
time <= end_time &
EUKAP2_00000018 > 0 &
animal_id=="12/10") %>%
lm(time ~ log(EUKAP2_00000018),data=.) -> lm_12_10
data_euka01 %>%
filter(time >= start_time &
time <= end_time &
EUKAP2_00000018 > 0) %>%
lm(time ~ log(EUKAP2_00000018),data=.) -> lm_all
half_life <- tibble(animal_id = c("9/10","10/10","12/10","all"),
licken_euka03 = -c(lm_9_10$coefficients[2],
lm_10_10$coefficients[2],
lm_12_10$coefficients[2],
lm_all$coefficients[2]) * log(2) * 24) %>%
mutate(ci_licken_euka03 = qnorm(p = c(0.975),
mean = 0,
sd = c(summary(lm_9_10)$coefficients[,"Std. Error"][2],
summary(lm_10_10)$coefficients[,"Std. Error"][2],
summary(lm_12_10)$coefficients[,"Std. Error"][2],
summary(lm_all)$coefficients[,"Std. Error"][2]))
* log(2) * 24,
ci_licken_euka03_inf = licken_euka03 - ci_licken_euka03,
ci_licken_euka03_sup = licken_euka03 + ci_licken_euka03)
half_lifeBetulaceae.taxid = ecofind(taxo,"^Betulaceae$")
is.a.Betulaceae = is.subcladeof(taxonomy = taxo,Euka01@motus$taxid,Betulaceae.taxid)
is.a.Betulaceae[is.na(is.a.Betulaceae)]=FALSE
sum(is.a.Betulaceae,na.rm = TRUE)## [1] 1Euka01@motus[is.a.Betulaceae,c("scientific_name","mean_ref_freq")]Betulaceae = Euka01@reads[,is.a.Lecanoromycetes]start_time=1
end_time=9
plot1 = ggplot(data = data_euka01,
aes(x = time,
y = EUKAP1_00000028,
color = animal_id)) +
geom_point() +
geom_smooth(data = data_euka01 %>% filter(time >= start_time & time <= end_time),
method = lm,show.legend = FALSE) +
scale_color_manual(values=cbbPalette) +
stat_summary_bin(fun = mean, geom = "line") +
scale_y_log10() +
geom_vline (xintercept = c(0.54,1.54,2.54), colour = cbbPalette[6]) +
geom_vline (xintercept = c(5.54,6.54,7.54), colour = cbbPalette[7]) +
geom_vline (xintercept = c(10.54,11.54,12.54), colour = cbbPalette[8]) +
geom_vline (xintercept = 26, colour = "black") +
scale_x_continuous(breaks = scales::pretty_breaks(n = 15),limits = c(-2,30)) +
guides(color=guide_legend(title="Individual")) +
theme_minimal() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black")) +
labs(y="Proportion of Betula DNA",
x="Time (days after Betula pubescens was fed)") +
labs(x=expression('Time (days after'~italic(Betula)~ ~italic( pubescens)~ 'was fed)'),
fill="Var1")
plot1## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 26 rows containing non-finite values (stat_summary_bin).## Warning: Removed 15 rows containing missing values (geom_point).
ggsave("Euka.Biomass_for_publication_lines.png", plot = plot1, width = 25, height = 16, units = c("cm"))## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 26 rows containing non-finite values (stat_summary_bin).## Warning: Removed 15 rows containing missing values (geom_point).data_euka01 %>%
filter(time >= start_time &
time <= end_time &
EUKAP1_00000028 > 0 &
animal_id=="9/10") %>%
lm(time ~ log(EUKAP1_00000028),data=.) -> lm_9_10
data_euka01 %>%
filter(time >= start_time &
time <= end_time &
EUKAP1_00000028 > 0 &
animal_id=="10/10") %>%
lm(time ~ log(EUKAP1_00000028),data=.) -> lm_10_10
data_euka01 %>%
filter(time >= start_time &
time <= end_time &
EUKAP1_00000028 > 0 &
animal_id=="12/10") %>%
lm(time ~ log(EUKAP1_00000028),data=.) -> lm_12_10
data_euka01 %>%
filter(time >= start_time &
time <= end_time &
EUKAP1_00000028 > 0) %>%
lm(time ~ log(EUKAP1_00000028),data=.) -> lm_all
half_life <- mutate(half_life,
animal_id = c("9/10","10/10","12/10","all"),
betula_euka03 = -c(lm_9_10$coefficients[2],
lm_10_10$coefficients[2],
lm_12_10$coefficients[2],
lm_all$coefficients[2]) * log(2) * 24,
ci_betula_euka03 = qnorm(p = c(0.975),
mean = 0,
sd = c(summary(lm_9_10)$coefficients[,"Std. Error"][2],
summary(lm_10_10)$coefficients[,"Std. Error"][2],
summary(lm_12_10)$coefficients[,"Std. Error"][2],
summary(lm_all)$coefficients[,"Std. Error"][2]))
* log(2) * 24,
ci_betula_euka03_inf = betula_euka03 - ci_betula_euka03,
ci_betula_euka03_sup = betula_euka03 + ci_betula_euka03)
half_lifeplot3 = ggplot(data = as.data.frame(Euka01@reads),
aes(x = samples(Euka01)$times_from_birch/24,
y = EUKAP2_00000018,
color = Animal)) +
geom_point(aes(color = Animal)) +
scale_x_continuous(breaks = c(0, 10, 20, 30)) +
xlim(-2,30) +
geom_vline (xintercept = 26, colour = "black") +
theme_minimal() +
theme(panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black")) +
labs(y="Proportion of Lecanoromycetidae DNA",
x="Time (days after Betula pubescens was fed)") +
labs(x=expression('Time (days after'~italic(Betula)~ ~italic( pubescens)~ 'was fed)'),
fill="Var1")## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.plot3## Warning: Removed 15 rows containing missing values (geom_point).
ggsave("Euka.Biomass_for_publication_scattercolor.png", plot = plot3, width = 25, height = 16, units = c("cm"))## Warning: Removed 15 rows containing missing values (geom_point).plot3 = ggplot(data = as.data.frame(Euka01@reads),
aes(x = samples(Euka01)$times_from_birch/24,
y = EUKAP2_00000018)) +
geom_point(aes(shape = Animal), size = 2) +
scale_shape_manual(values=c(4, 2, 1))+
scale_x_continuous(breaks = c(0, 10, 20, 30)) +
xlim(-2,30) +
geom_vline (xintercept = 26, colour = "black") +
theme_minimal() +
theme(panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black")) +
labs(y="Proportion of Lecanoromycetidae DNA",
x="Time (days after Betula pubescens was fed)") +
labs(x=expression('Time (days after'~italic(Betula)~ ~italic( pubescens)~ 'was fed)'),
fill="Var1")## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.plot3## Warning: Removed 15 rows containing missing values (geom_point).
plot1 = ggplot(data = as.data.frame(Euka01@reads),
aes(x = samples(Euka01)$times_from_birch/24,
y = EUKAP2_00000018,
color = Animal)) +
scale_x_continuous(breaks = c(0, 10, 20, 30)) +
stat_summary_bin(fun.y = mean, geom = "line") +
xlim(-2,30) +
geom_vline (xintercept = 0.54, colour = "lightgrey") +
geom_vline (xintercept = 1.54, colour = "lightgrey") +
geom_vline (xintercept = 2.54, colour = "lightgrey") +
geom_vline (xintercept = 5.54, colour = "lightgrey") +
geom_vline (xintercept = 6.54, colour = "lightgrey") +
geom_vline (xintercept = 7.54, colour = "lightgrey") +
geom_vline (xintercept = 10.54, colour = "lightgrey") +
geom_vline (xintercept = 11.54, colour = "lightgrey") +
geom_vline (xintercept = 12.54, colour = "lightgrey") +
annotate("text", x = 1.60, y = 0.23, label = "20g") +
annotate("text", x = 6.60, y = 0.23, label = "500g") +
annotate("text", x = 11.60, y = 0.23, label = "2000g") +
geom_vline (xintercept = 26, colour = "black") +
theme_minimal() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black")) +
labs(y="Proportion of Lecanoromycetidae DNA",
x="Time (days after Betula pubescens was fed)") +
labs(x=expression('Time (days after'~italic(Betula)~ ~italic( pubescens)~ 'was fed)'),
fill="Var1")## Warning: `fun.y` is deprecated. Use `fun` instead.## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.plot1## Warning: Removed 15 rows containing non-finite values (stat_summary_bin).
plot1 = ggplot(data = as.data.frame(Euka01@reads),
aes(x = samples(Euka01)$times_from_birch/24,
y = EUKAP2_00000018)) +
stat_summary_bin(fun.y = mean, geom = "line") +
geom_point(aes(color = Animal)) +
xlim(-2,30) +
geom_vline (xintercept = 0.54, colour = "lightgrey") +
geom_vline (xintercept = 1.54, colour = "lightgrey") +
geom_vline (xintercept = 2.54, colour = "lightgrey") +
geom_vline (xintercept = 5.54, colour = "lightgrey") +
geom_vline (xintercept = 6.54, colour = "lightgrey") +
geom_vline (xintercept = 7.54, colour = "lightgrey") +
geom_vline (xintercept = 10.54, colour = "lightgrey") +
geom_vline (xintercept = 11.54, colour = "lightgrey") +
geom_vline (xintercept = 12.54, colour = "lightgrey") +
annotate("text", x = 1.60, y = 0.27, label = "20g") +
annotate("text", x = 6.60, y = 0.27, label = "500g") +
annotate("text", x = 11.60, y = 0.27, label = "2000g") +
geom_vline (xintercept = 25, colour = "black", linetype = "dotted") +
theme_minimal() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black")) +
labs(y="Proportion of lichen DNA",
x="Time (days after B. pubescens was fed)") +
labs(x=expression('Time (days after'~italic(B.)~ ~italic( pubescens)~ 'was fed)'),
fill="Var1")## Warning: `fun.y` is deprecated. Use `fun` instead.plot1## Warning: Removed 15 rows containing non-finite values (stat_summary_bin).## Warning: Removed 15 rows containing missing values (geom_point).
ggsave("Euka.Biomass_for_publication_onelineanddots.png", plot = plot1, width = 25, height = 16, units = c("cm"))## Warning: Removed 15 rows containing non-finite values (stat_summary_bin).## Warning: Removed 15 rows containing missing values (geom_point).data = cbind(as.data.frame(Euka01@reads),
time = samples(Euka01)$times_from_birch/24,
animal_id = samples(Euka01)$animal_id,
amount = 2000) %>%
mutate(amount = ifelse(time < 10,500,amount)) %>%
mutate(amount = ifelse(time < 5,20,amount))
start_time=11
end_time=23
plot1 = ggplot(data = data,
aes(x = time,
y = EUKAP2_00000018,
col=animal_id)) +
geom_smooth(span=0.2) +
geom_point() +
scale_y_log10() +
geom_vline (xintercept = c(0.54,1.54,2.54), colour = cbbPalette[6]) +
geom_vline (xintercept = c(5.54,6.54,7.54), colour = cbbPalette[7]) +
geom_vline (xintercept = c(10.54,11.54,12.54), colour = cbbPalette[8]) +
geom_vline (xintercept = 26, colour = "black") +
guides(color=guide_legend(title="Individual")) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 15),limits = c(-2,30)) +
theme_minimal() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black")) +
labs(y="Proportion of Lecanoromycetidae DNA",
x="Time (days after Betula pubescens was fed)") +
labs(x=expression('Time (days after'~italic(Betula)~ ~italic( pubescens)~ 'was fed)'),
fill="Var1")
plot1
data %>% group_by(amount,animal_id) %>% summarise(max_lichen = max(EUKAP2_00000018)) %>%
ggplot(aes(x=amount,y=max_lichen)) + geom_point(aes(col=animal_id))
data %>%
group_by(animal_id) %>%
mutate(lichen=EUKAP2_00000018,
sm = rollapply(EUKAP2_00000018,3,
mean,align="center",
fill=NA)) %>%
select(-starts_with("EUKA")) %>%
group_by(animal_id,amount) %>%
summarise(max_lichen = max(sm,na.rm = TRUE),
m2 = max(lichen)) %>%
ggplot(aes(x=amount,y=m2)) +
geom_point(aes(col=animal_id)) +
scale_x_log10() +
scale_y_log10() +
geom_smooth(method=lm)## `summarise()` has grouped output by 'animal_id'. You can override using the `.groups` argument.## `geom_smooth()` using formula 'y ~ x'
data %>%
group_by(animal_id) %>%
mutate(lichen=EUKAP2_00000018,
sm = rollapply(EUKAP2_00000018,3,
mean,align="center",
fill=NA)) %>%
select(-starts_with("EUKA")) %>%
group_by(animal_id,amount) %>%
summarise(max_lichen = max(sm,na.rm = TRUE),
m2 = max(lichen),
time = rev(time[which.max(lichen)])[1]) %>%
mutate(last_feed = ifelse(amount==20,2.54,ifelse(amount==500,7.54,12.54)),
delta = time - last_feed)## `summarise()` has grouped output by 'animal_id'. You can override using the `.groups` argument.data %>% group_by(animal_id) %>% mutate(lichen=EUKAP2_00000018,sm = rollapply(EUKAP2_00000018,3,mean,align="center",fill=NA)) %>% select(-starts_with("EUKA")) %>% group_by(animal_id,amount) %>% summarise(max_lichen = max(sm,na.rm = TRUE),m2 = max(lichen)) %>%
lm(log(amount) ~ log(m2),data = .) -> lm_gram_rra## `summarise()` has grouped output by 'animal_id'. You can override using the `.groups` argument.summary(lm_gram_rra)##
## Call:
## lm(formula = log(amount) ~ log(m2), data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.08798 -0.00749 -0.00098 0.16311 0.61214
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.8327 0.4354 22.58 8.45e-08 ***
## log(m2) 1.5378 0.1448 10.62 1.43e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5287 on 7 degrees of freedom
## Multiple R-squared: 0.9416, Adjusted R-squared: 0.9332
## F-statistic: 112.8 on 1 and 7 DF, p-value: 1.435e-05data$lichen_gr <- data$EUKAP2_00000018^coefficients(lm_gram_rra)[2] * exp(coefficients(lm_gram_rra)[1])data %>% ggplot(aes(x = lichen_gr,
y = EUKAP2_00000018,
col=animal_id)) + geom_point()
start_time=11
end_time=23
plot1 = ggplot(data = data,
aes(x = time,
y = lichen_gr,
col=animal_id)) +
geom_smooth(span=0.2) +
geom_point() +
scale_y_log10() +
geom_vline (xintercept = c(0.54,1.54,2.54), colour = cbbPalette[6]) +
geom_vline (xintercept = c(5.54,6.54,7.54), colour = cbbPalette[7]) +
geom_vline (xintercept = c(10.54,11.54,12.54), colour = cbbPalette[8]) +
geom_vline (xintercept = 26, colour = "black") +
guides(color=guide_legend(title="Individual")) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 15),limits = c(-2,30)) +
theme_minimal() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black")) +
labs(y="Lichen (g)",
x="Time (days after Betula pubescens was fed)") +
labs(x=expression('Time (days after'~italic(Betula)~ ~italic( pubescens)~ 'was fed)'),
fill="Var1")
plot1
data %>% group_by(amount,animal_id) %>% summarise(max_lichen = max(EUKAP2_00000018)) %>%
ggplot(aes(x=amount,y=max_lichen)) + geom_point(aes(col=animal_id))
data %>%
filter(time >= start_time &
time <= end_time &
lichen_gr > 0 &
animal_id=="9/10") %>%
lm(time ~ log(lichen_gr),data=.) -> lm_9_10
data %>%
filter(time >= start_time &
time <= end_time &
lichen_gr > 0 &
animal_id=="10/10") %>%
lm(time ~ log(lichen_gr),data=.) -> lm_10_10
data %>%
filter(time >= start_time &
time <= end_time &
lichen_gr > 0 &
animal_id=="12/10") %>%
lm(time ~ log(lichen_gr),data=.) -> lm_12_10
data %>%
filter(time >= start_time &
time <= end_time &
lichen_gr > 0) %>%
lm(time ~ log(lichen_gr),data=.) -> lm_all
half_life <- tibble(animal_id = c("9/10","10/10","12/10","all"),
licken_euka03 = -c(lm_9_10$coefficients[2],
lm_10_10$coefficients[2],
lm_12_10$coefficients[2],
lm_all$coefficients[2]) * log(2) * 24) %>%
mutate(ci_licken_euka03 = qnorm(p = c(0.975),
mean = 0,
sd = c(summary(lm_9_10)$coefficients[,"Std. Error"][2],
summary(lm_10_10)$coefficients[,"Std. Error"][2],
summary(lm_12_10)$coefficients[,"Std. Error"][2],
summary(lm_all)$coefficients[,"Std. Error"][2]))
* log(2) * 24,
ci_licken_euka03_inf = licken_euka03 - ci_licken_euka03,
ci_licken_euka03_sup = licken_euka03 + ci_licken_euka03)
half_life\[ log(RRA) = a \times log(lichen) + b \] \[ RRA = e^b \times lichen^a \]
Euka01_plant <- Euka01[,Euka01@motus$is_plant]
dim(Euka01_plant)## [1] 210 17Euka01_plant@reads %>%
as_tibble() %>%
mutate(sample_id = rownames(Euka01_plant)) %>%
pivot_longer(cols = - "sample_id", names_to = "motu_id",values_to = "relfreq") %>%
left_join(Euka01_plant@motus, by = c(motu_id = "id")) %>%
group_by(sample_id,family_name) %>%
summarise(relfreq = sum(relfreq)) %>%
ggplot(aes(x=family_name,y=relfreq)) +
geom_boxplot() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))## `summarise()` has grouped output by 'sample_id'. You can override using the `.groups` argument.
Euka01_plant@reads %>%
as_tibble() %>%
mutate(sample_id = rownames(Euka01_plant)) %>%
pivot_longer(cols = - "sample_id", names_to = "motu_id",values_to = "relfreq") %>%
left_join(Euka01_plant@motus, by = c(motu_id = "id")) %>%
group_by(sample_id,family_name) %>%
summarise(relfreq = sum(relfreq)) %>%
ungroup() %>%
left_join(Euka01_plant@samples, by = "sample_id") %>%
ggplot(aes(x=times_from_birch/24,y=relfreq,col=family_name)) +
geom_point( ) +
xlim(-2,30) +
facet_grid(family_name ~ .,scales="free_y") +
geom_vline (xintercept = 0.54, colour = "lightgrey") +
geom_vline (xintercept = 1.54, colour = "lightgrey") +
geom_vline (xintercept = 2.54, colour = "lightgrey") +
geom_vline (xintercept = 5.54, colour = "lightgrey") +
geom_vline (xintercept = 6.54, colour = "lightgrey") +
geom_vline (xintercept = 7.54, colour = "lightgrey") +
geom_vline (xintercept = 10.54, colour = "lightgrey") +
geom_vline (xintercept = 11.54, colour = "lightgrey") +
geom_vline (xintercept = 12.54, colour = "lightgrey") +
annotate("text", x = 1.60, y = 0.27, label = "20g") +
annotate("text", x = 6.60, y = 0.27, label = "500g") +
annotate("text", x = 11.60, y = 0.27, label = "2000g") +
geom_vline (xintercept = 25, colour = "black", linetype = "dotted") ## `summarise()` has grouped output by 'sample_id'. You can override using the `.groups` argument.## Warning: Removed 135 rows containing missing values (geom_point).
Attempt to correct relative frequencies by the average relative frequency of constantes plants
Euka01_plant@reads %>%
as_tibble() %>%
mutate(sample_id = rownames(Euka01_plant)) %>%
pivot_longer(cols = - "sample_id", names_to = "motu_id",values_to = "relfreq") %>%
left_join(Euka01_plant@motus %>% select(- starts_with("obiclean_status") ),
by = c(motu_id = "id")) %>%
ungroup() %>%
left_join(Euka01_plant@samples, by = "sample_id") %>%
mutate(times_from_birch = times_from_birch/24,
time_group = floor(times_from_birch)) %>%
filter(time_group < 26) %>%
group_by(family_name) %>%
summarise(mean_relfreq = mean(relfreq)) %>%
arrange(desc(mean_relfreq)) -> family_means
family_meansEuka01_plant@reads %>%
as_tibble() %>%
mutate(sample_id = rownames(Euka01_plant)) %>%
pivot_longer(cols = - "sample_id", names_to = "motu_id",values_to = "relfreq") %>%
left_join(Euka01_plant@motus %>% select(- starts_with("obiclean_status") ),
by = c(motu_id = "id")) %>%
ungroup() %>%
left_join(Euka01_plant@samples, by = "sample_id") %>%
mutate(times_from_birch = times_from_birch/24,
time_group = floor(times_from_birch)) %>%
filter(time_group < 26) %>%
group_by(family_name,animal_id) %>%
mutate(mean_relfreq = mean(relfreq)) %>%
arrange(desc(mean_relfreq)) %>%
group_by(family_name,animal_id,mean_relfreq,time_group) %>%
summarise(mean_group = mean(relfreq)) %>%
mutate(correction = mean_relfreq/mean_group) %>%
ungroup() -> corrections## `summarise()` has grouped output by 'family_name', 'animal_id', 'mean_relfreq'. You can override using the `.groups` argument.correctionsstart_time=1
end_time=9
Euka01@reads %>%
as_tibble() %>%
mutate(sample_id = rownames(Euka01)) %>%
pivot_longer(cols = - "sample_id", names_to = "motu_id",values_to = "relfreq") %>%
left_join(Euka01@motus %>% select(- starts_with("obiclean_status") ),
by = c(motu_id = "id")) %>%
ungroup() %>%
left_join(Euka01@samples, by = "sample_id") %>%
mutate(times_from_birch = times_from_birch/24,
time_group = floor(times_from_birch)) %>%
filter(time_group < 26) %>%
left_join(corrections %>%
filter(family_name == "Brassicaceae") %>%
select(Animal_id = animal_id, time_group,correction)) %>%
mutate(cor_relfreq = relfreq * correction) %>%
filter(family_name=="Betulaceae") %>%
group_by(sample_id,Animal_id,times_from_birch,time_group) %>%
summarise(relfreq=sum(cor_relfreq)) -> corrrected_data_Euka01## Joining, by = "time_group"## `summarise()` has grouped output by 'sample_id', 'Animal_id', 'times_from_birch'. You can override using the `.groups` argument.ggplot(data = corrrected_data_Euka01,
aes(x = times_from_birch,
y = relfreq,
color = Animal_id)) +
geom_point() +
geom_smooth(data = corrrected_data_Euka01 %>% filter(times_from_birch >= start_time & times_from_birch <= end_time),
method = lm,show.legend = FALSE) +
scale_color_manual(values=cbbPalette) +
stat_summary_bin(fun = mean, geom = "line") +
scale_y_log10() +
geom_vline (xintercept = c(0.54,1.54,2.54), colour = cbbPalette[6]) +
geom_vline (xintercept = c(5.54,6.54,7.54), colour = cbbPalette[7]) +
geom_vline (xintercept = c(10.54,11.54,12.54), colour = cbbPalette[8]) +
geom_vline (xintercept = 26, colour = "black") +
scale_x_continuous(breaks = scales::pretty_breaks(n = 15),limits = c(-2,30)) +
guides(color=guide_legend(title="Individual")) +
theme_minimal() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black")) +
labs(y="Proportion of Betula DNA",
x="Time (days after Betula pubescens was fed)") +
labs(x=expression('Time (days after'~italic(Betula)~ ~italic( pubescens)~ 'was fed)'),
fill="Var1")## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 33 rows containing non-finite values (stat_summary_bin).## Warning: Removed 1 rows containing missing values (geom_point).
corrrected_data_Euka01 %>%
filter(times_from_birch >= start_time &
times_from_birch <= end_time &
relfreq > 0 &
Animal_id=="9/10") %>%
lm(times_from_birch ~ log(relfreq),data=.) -> lm_9_10
corrrected_data_Euka01 %>%
filter(times_from_birch >= start_time &
times_from_birch <= end_time &
relfreq > 0 &
Animal_id=="10/10") %>%
lm(times_from_birch ~ log(relfreq),data=.) -> lm_10_10
corrrected_data_Euka01 %>%
filter(times_from_birch >= start_time &
times_from_birch <= end_time &
relfreq > 0 &
Animal_id=="12/10") %>%
lm(times_from_birch ~ log(relfreq),data=.) -> lm_12_10
corrrected_data_Euka01 %>%
filter(times_from_birch >= start_time &
times_from_birch <= end_time &
relfreq > 0) %>%
lm(times_from_birch ~ log(relfreq) + Animal_id,data=.) -> lm_all
half_life <- tibble(animal_id = c("9/10","10/10","12/10","all"),
betula_euka02 = -c(lm_9_10$coefficients[2],
lm_10_10$coefficients[2],
lm_12_10$coefficients[2],
lm_all$coefficients[2]) * log(2) * 24,
ci_betula_euka02 = qnorm(p = c(0.975),
mean = 0,
sd = c(summary(lm_9_10)$coefficients[,"Std. Error"][2],
summary(lm_10_10)$coefficients[,"Std. Error"][2],
summary(lm_12_10)$coefficients[,"Std. Error"][2],
summary(lm_all)$coefficients[,"Std. Error"][2]))
* log(2) * 24,
ci_betula_euka02_inf = betula_euka02 - ci_betula_euka02,
ci_betula_euka02_sup = betula_euka02 + ci_betula_euka02)
half_lifeFor lickens
start_time=11
end_time=23
Euka01@reads %>%
as_tibble() %>%
mutate(sample_id = rownames(Euka01)) %>%
pivot_longer(cols = - "sample_id", names_to = "motu_id",values_to = "relfreq") %>%
left_join(Euka01@motus %>% select(- starts_with("obiclean_status") ),
by = c(motu_id = "id")) %>%
ungroup() %>%
left_join(Euka01@samples, by = "sample_id") %>%
mutate(times_from_birch = times_from_birch/24,
time_group = floor(times_from_birch)) %>%
filter(time_group < 26) %>%
left_join(corrections %>%
filter(family_name == "Brassicaceae") %>%
select(Animal_id = animal_id, time_group,correction)) %>%
mutate(cor_relfreq = relfreq * correction) %>%
filter(motu_id=="EUKAP2_00000018") %>%
group_by(sample_id,Animal_id,times_from_birch,time_group) %>%
summarise(relfreq=sum(cor_relfreq)) -> corrrected_data_Euka01## Joining, by = "time_group"## `summarise()` has grouped output by 'sample_id', 'Animal_id', 'times_from_birch'. You can override using the `.groups` argument.ggplot(data = corrrected_data_Euka01,
aes(x = times_from_birch,
y = relfreq,
color = Animal_id)) +
geom_point() +
geom_smooth(data = corrrected_data_Euka01 %>% filter(times_from_birch >= start_time & times_from_birch <= end_time),
method = lm,show.legend = FALSE) +
scale_color_manual(values=cbbPalette) +
stat_summary_bin(fun = mean, geom = "line") +
scale_y_log10() +
geom_vline (xintercept = c(0.54,1.54,2.54), colour = cbbPalette[6]) +
geom_vline (xintercept = c(5.54,6.54,7.54), colour = cbbPalette[7]) +
geom_vline (xintercept = c(10.54,11.54,12.54), colour = cbbPalette[8]) +
geom_vline (xintercept = 26, colour = "black") +
scale_x_continuous(breaks = scales::pretty_breaks(n = 15),limits = c(-2,30)) +
guides(color=guide_legend(title="Individual")) +
theme_minimal() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black")) +
labs(y="Proportion of Betula DNA",
x="Time (days after Betula pubescens was fed)") +
labs(x=expression('Time (days after'~italic(Betula)~ ~italic( pubescens)~ 'was fed)'),
fill="Var1")## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 24 rows containing non-finite values (stat_smooth).## Warning: Removed 66 rows containing non-finite values (stat_summary_bin).## Warning: Removed 1 rows containing missing values (geom_point).
corrrected_data_Euka01 %>%
filter(times_from_birch >= start_time &
times_from_birch <= end_time &
relfreq > 0 &
Animal_id=="9/10") %>%
lm(times_from_birch ~ log(relfreq),data=.) -> lm_9_10
corrrected_data_Euka01 %>%
filter(times_from_birch >= start_time &
times_from_birch <= end_time &
relfreq > 0 &
Animal_id=="10/10") %>%
lm(times_from_birch ~ log(relfreq),data=.) -> lm_10_10
corrrected_data_Euka01 %>%
filter(times_from_birch >= start_time &
times_from_birch <= end_time &
relfreq > 0 & relfreq <= 1,
Animal_id=="12/10") %>%
lm(times_from_birch ~ log(relfreq),data=.) -> lm_12_10
corrrected_data_Euka01 %>%
filter(times_from_birch >= start_time &
times_from_birch <= end_time &
relfreq > 0& relfreq <= 1) %>%
lm(times_from_birch ~ log(relfreq) + Animal_id,data=.) -> lm_all
half_life <- tibble(animal_id = c("9/10","10/10","12/10","all"),
lichen_euka02 = -c(lm_9_10$coefficients[2],
lm_10_10$coefficients[2],
lm_12_10$coefficients[2],
lm_all$coefficients[2]) * log(2) * 24,
ci_lichen_euka02 = qnorm(p = c(0.975),
mean = 0,
sd = c(summary(lm_9_10)$coefficients[,"Std. Error"][2],
summary(lm_10_10)$coefficients[,"Std. Error"][2],
summary(lm_12_10)$coefficients[,"Std. Error"][2],
summary(lm_all)$coefficients[,"Std. Error"][2]))
* log(2) * 24,
ci_lichen_euka02_inf = lichen_euka02 - ci_lichen_euka02,
ci_lichen_euka02_sup = lichen_euka02 + ci_lichen_euka02)
half_lifeEffect on the relatiionship with the amount of ingested lichens
data = corrrected_data_Euka01 %>%
mutate(time = times_from_birch,
amount = 2000) %>%
mutate(amount = ifelse(time < 10,500,amount)) %>%
mutate(amount = ifelse(time < 5,20,amount))
start_time=11
end_time=23
plot1 = ggplot(data = data,
aes(x = time,
y = relfreq,
col=Animal_id)) +
geom_smooth(span=0.2) +
geom_point() +
scale_y_log10() +
geom_vline (xintercept = c(0.54,1.54,2.54), colour = cbbPalette[6]) +
geom_vline (xintercept = c(5.54,6.54,7.54), colour = cbbPalette[7]) +
geom_vline (xintercept = c(10.54,11.54,12.54), colour = cbbPalette[8]) +
geom_vline (xintercept = 26, colour = "black") +
guides(color=guide_legend(title="Individual")) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 15),limits = c(-2,30)) +
theme_minimal() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black")) +
labs(y="Proportion of Lecanoromycetidae DNA",
x="Time (days after Betula pubescens was fed)") +
labs(x=expression('Time (days after'~italic(Betula)~ ~italic( pubescens)~ 'was fed)'),
fill="Var1")
plot1
data %>% group_by(amount,Animal_id) %>% summarise(max_lichen = max(relfreq, na.rm = TRUE)) %>%
ggplot(aes(x=amount,y=max_lichen)) + geom_point(aes(col=Animal_id))
data %>%
group_by(Animal_id) %>%
mutate(lichen=relfreq) %>%
group_by(Animal_id,amount) %>%
summarise(max_lichen = max(lichen)) %>%
ggplot(aes(x=amount,y=max_lichen)) +
geom_point(aes(col=Animal_id)) +
scale_x_log10() +
scale_y_log10() +
geom_smooth(method=lm)## `summarise()` has grouped output by 'Animal_id'. You can override using the `.groups` argument.## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 1 rows containing non-finite values (stat_smooth).## Warning: Removed 1 rows containing missing values (geom_point).
data %>% group_by(Animal_id) %>%
mutate(lichen=relfreq) %>%
group_by(Animal_id,amount) %>%
summarise(max_lichen = max(lichen)) %>%
lm(log(amount) ~ log(max_lichen),data = .) -> lm_gram_rra## `summarise()` has grouped output by 'Animal_id'. You can override using the `.groups` argument.summary(lm_gram_rra)##
## Call:
## lm(formula = log(amount) ~ log(max_lichen), data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.96132 -0.55852 -0.07715 0.51804 1.15159
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.3342 0.4286 17.111 2.55e-06 ***
## log(max_lichen) 1.0419 0.1681 6.198 0.000813 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8081 on 6 degrees of freedom
## (1 observation effacée parce que manquante)
## Multiple R-squared: 0.8649, Adjusted R-squared: 0.8424
## F-statistic: 38.41 on 1 and 6 DF, p-value: 0.0008132